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Quadrature  moments  method  for  the  simulation 
of  turbulent  reactive  flows 

By  Venkatramanan  Raman,  Heinz  Pitsch  and  Rodney  O.  Fox  f 


1.  Motivation  and  Objectives 

In  recent  years,  computational  fluid  mechanics  has  become  one  of  the  primary  tools  for 
design  and  optimization  of  chemical  reactors.  With  stringent  environmental  constraints, 
close  control  of  product  selectivity  and  an  estimate  of  by-products  are  essential  in  suc- 
cessfully operating  chemical  plants.  The  fast  throughput  and  enhanced  mixing  conditions 
offered  by  turbulent  flows  are  increasingly  exploited  in  chemical  reactors.  Viable  simula- 
tion methods  for  such  flows  should  be  able  to  model  the  complex  interaction  of  reaction 
and  turbulent  flow.  All  known  reaction  models  used  in  these  simulations  follow  a segre- 
gated approach  where  different  techniques  are  used  to  solve  the  momentum  and  scalar 
transport  equations.  It  is  assumed  that  reaction  affects  fluid  flow  only  through  the  change 
in  density.  Usually  a variable  density  flow  solver  is  utilized  that  accepts  the  density  field 
from  the  scalar  handler  to  correct  the  flow  field.  The  scalar  transport  scheme  uses  the  flow 
properties  to  evaluate  the  local  density  and  this  iterative  procedure  is  used  to  advance 
the  solution  in  time.  The  Eulerian  solution  technique  that  is  commonly  used  in  solving 
scalar  transport  equations  inherently  does  not  contain  information  about  sub-grid  level 
processes.  The  adverse  effect  of  neglecting  sub-grid  scalar  fluctuations  in  cases  where 
the  scalar  evolves  through  non-linear  rate  expressions  is  well  known.  In  combustion  pro- 
cesses which  are  characterized  by  fast  chemistry,  use  of  fiamelet  model  (Pitsch  & Steiner 
2000)  or  conditional  moment  closure  (CMC)  (Bilger  1993)  based  on  a conserved  scalar 
is  known  to  be  quite  accurate  in  making  qualitative  as  well  as  quantitative  predictions. 
Such  a method  obviates  the  need  for  solving  multiple  scalar  transport  equations  and 
is  not  restricted  by  the  time-scale  of  individual  reactions  in  the  chemistry  mechanism. 
The  fiamelet  model,  like  all  other  reaction  models,  requires  the  specification  of  a scalar 
dissipation  rate  and  assumes  the  shape  of  the  PDF  at  the  sub-grid  level.  The  first-order 
CMC  method  models  the  conditional  mean  of  the  reacting  scalars  as  a function  of  mix- 
ture fraction  and  local  flow  conditions.  But  the  model  ignores  any  fluctuations  about  the 
conditional  mean  and  may  not  be  viable  in  slow  chemistry  regimes.  Both  the  CMC  and 
fiamelet  model  assume  that  all  the  reactive  scalars  can  be  parameterized  by  a single  con- 
served scalar.  On  the  other  end  of  the  spectrum,  the  transported-PDF  scheme  computes 
the  sub-grid  scalar-PDF  in  terms  of  a set  of  delta-functions.  This  method  can  be  used 
in  tandem  with  a flow  solver  like  those  based  on  the  Reynolds- Averaged  Navier  Stokes 
(RANS)  equation  or  Large-Eddy  Simulation  (LES)  scheme  to  model  reaction  (Muradoglu 
et  al.  1999;  Haworth  & El  Tahry  1991).  Such  a formulation  computes  the  profiles  of  all 
the  species  involved  in  the  flow  and  leads  to  a closed  form  for  the  reaction  source  term. 
However,  a mixing  model  is  needed  to  describe  the  sub-grid  mixing  process  and  has  been 
the  focus  of  study  in  PDF  methods  (Subramaniam  & Pope  1999;  Valino  1995;  Tsai  & Fox 
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1998).  The  fact  that  no  convincing  mixing  model  exists  is  one  of  the  important  drawbacks 
of  the  method. 

The  transported-PDF  still  is  an  attractive  scheme  in  that  the  reaction  source  term 
appears  closed.  The  inhomogeneous  transport  equation  for  the  sub-grid  PDF  is  multi- 
dimensional and  cannot  be  solved  using  Eulerian  grid  techniques.  Particle  based  Monte- 
Carlo  schemes  used  with  realistic  chemistry  show  good  agreement  with  experimental 
results  (Xu  & Pope  2000;  Masri  & Pope  1990;  Roekaerts  1991;  Raman  et  al.  2003a). 
However,  the  handling  of  a large  number  of  Lagrangian  particles  along  with  detailed 
chemistry  can  be  computationally  expensive.  When  using  LES  as  the  flow  solver,  such 
schemes  can  be  outright  intractable  in  practical  flow  configurations.  In  this  work,  an 
equivalent  Eulerian  version  of  the  transported-PDF  method  is  formulated  for  LES  of 
reactive  flows. 


2.  Filtered  Quadrature  Scheme  for  Variable  Density  Flows 

For  use  in  LES  of  reactive  flows,  a filtered  density  function  (FDF)  is  defined  (Gao 
& O’Brien  1993)  that  prescribes  the  sub-filter  joint  composition  density  function  of  the 
scalars. 

/+oo 

p(x\  t)5[ij>  - 0(x',  t))G(x'  - x)dx',  (2.1) 

-OO 

where  Fl  is  the  mass  density  function.  This  definition  implies  that  in  variable  density 
flows,  the  FDF  is  the  mass  weighted  and  spatially  filtered  fine-grain  density  (Jaberi  et  al. 
1999). 
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FL(i>-,x,t)dxp  = p,  (2.2) 
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where  p is  the  filtered  density.  Using  the  scalar  transport  equation,  the  FDF  is  shown  to 
evolve  in  multi-dimensional  space  as: 
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where  is  the  Favre  filtered  velocity  component,  ip  is  the  composition  array  and  (p 
is  the  Favre  filtered  mean  composition  array.  The  molecular  and  turbulent  diffusivi- 
ties  are  denoted  by  7 and  7 1 respectively.  The  mixing  frequency  is  modeled  as  flm  = 
Or  (7  + 7t)  / ({p)l  Aq),  where  Ac  is  the  filter  width.  This  high-  dimensional  equation 
is  traditionally  solved  using  Monte-Carlo  schemes  (Colucci  et  al.  1998;  Raman  et  al. 
20036).  As  mentioned  before,  the  particle  based  Monte-Carlo  schemes  are  computation- 
ally expensive  for  large  grids  and  hence  become  intractable  even  for  simple  chemistry. 
The  Direct  Quadrature  Method  of  Moments  (DQMOM)  (Marchisio  & Fox  2003;  Wang 
& Fox  2003)  is  introduced  here  in  the  context  of  LES  to  provide  an  alternate  tractable 
scheme. 

Figure  1 shows  the  fine-grained  PDF  composed  of  a set  of  delta-functions  approximat- 
ing an  arbitrary  shaped  PDF  at  a given  point  in  the  flow.  Any  particle  based  solution 
to  the  FDF  equation  will  yield  such  an  approximation.  The  FDF  plot  can  be  considered 
as  a plot  of  normalized  weights  of  particles  in  the  computational  cell.  The  accuracy  of 
the  approximation  depends,  among  other  factors,  on  the  number  of  approximating  delta 
functions.  The  Na  delta-functions  are  characterized  by  their  positions  (<pai)  and  their 
heights  ( u>ai ).  The  FDF  is  considered  solved  if  for  a given  number  of  delta  functions,  the 
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FIGURE  1 . Approximation  of  a PDF  using  a finite  number  of  delta  functions.  (Left) 
Transported  PDF  method  and  (Right)  DQMOM  method. 


values  of  <j>ai  and  wai  are  known.  The  DQMOM  model  formulates  transport  equations 
for  these  quantities  based  on  pre-determined  number  of  so-called  environments,  where 
each  environment  corresponds  to  a single  delta-peak.  The  source  terms  involved  in  these 
equations  ensure  that  the  DQMOM  equations  match  the  moment  equations  of  the  scalar. 


In  the  multi-environment  DQMOM  model,  the  FDF  is  assumed  to  be  of  the  form: 


N 

Fl  =p'%2wi8(il>-4>i).  (2.4) 

i=l 

The  individual  transport  equations  for  the  weights  and  locations  of  the  delta-functions 
can  be  derived  by  substituting  Eq.  2.4  into  Eq.  2.3  (Fox  2003).  For  a single  scalar  case, 
this  leads  to 
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where  G*  = and  R(ip\  x,  t)  contains  the  mixing  and  reaction  terms.  6^  indicates  the 
m-th  derivative  of  the  delta  function.  The  properties  of  5 function  are  used  in  rewriting 
the  derivatives  in  terms  of  the  filtered  quantities  (Pope  2000). 
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dx? 


(2.7) 


Following  Fox  (2003),  Eq.  2.5  can  be  rewritten  in  terms  of  the  transport  equations  for 
Wi  and  Gi,  using  a,  and  bi  to  denote  the  respective  equations. 


T:  [<$  (V> - d>)  + <M  1 (ip - d>)]  ai-^2 51  ~^)bi 

1=1  L 1=1 
N 

= s(2)  (ip-<p)  WiCi  + R(if>- x,  t), 


(2.8) 
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where  c<  = (r  + Ft)(Vdi)2.  In  deriving  these  equations,  the  only  assumption  made  so 
far  is  that  the  shape  of  the  PDF  is  approximated  by  a finite-set  of  delta  functions.  For 
a iV-environment  model,  2 N source  terms  need  to  be  specified.  It  should  be  noted  that 
if  the  number  of  scalars  is  more  than  one,  several  cross  moments  can  also  be  specified. 
For  details  on  the  feasibility  and  choice  of  moments  refer  to  Fox  (2003)  and  Marchisio 
&:  Fox  (2003).  Here,  only  the  pure  moments,  namely  the  mean  and  variance  of  a scalar 
are  used  in  fixing  the  source  terms.  Multiplying  Eq.  2.8  by  ipm  and  integrating  over  the 
composition  space  yields: 
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Similarly  the  reaction  and  mixing  terms  can  be  integrated  to  yield: 
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Using  Eqs  2.9  and  2.10  we  can  obtain  the  source  terms  ai  and  bi  by  solving  a set  of  non- 
linear algberaic  equations.  In  this  work,  N is  set  to  two  for  all  simulations.  To  further 
simplify  the  equations,  we  can  arbitrarily  set  a*  to  zero  (Fox  2003).  We  then  use  m — 1,2 
to  obtain  the  other  source  terms.  Using  the  above  moment  equations,  the  non-linear 
system  for  bi  reduces  to: 
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Prom  the  above  set  of  equations,  the  source  terms  are  determined  to  be: 

CL  1 ==  &2  ~ 0,  (2.12) 

1 2 

h = = ■=-  T WiTtiV^)2  + Om  hG2  - w2G1)  + w1S(<f>1),  (2.13) 

Vl  — V 2 t=l 

and 

i 2 

h = = =-  T WiTtiW^)2  - Clm  (w\G2  - w2G{)  + w2S(<t>2).  (2.14) 

01  ~ 02  i=j 

As  discussed  earlier,  the  number  of  environments  can  be  considered  to  be  equivalent  to 
the  particle  number  density  used  in  a Lagrangian  simulation.  This  implies  that  as  the 
number  of  environments  is  increased,  the  accuracy  of  the  scheme  should  increase  (Wang 
& Fox  2003).  The  single  scalar  DQMOM  can  be  extended  to  a multi- variate  case  using 
similar  derivation  technique  (Fox  2003).  In  determining  the  source  terms,  one  of  the 
important  criteria  is  that  the  pre-multiplier  of  the  b vector  in  Eq.  2.11  is  invertible.  It 
should  also  be  noted  that  a*  need  not  be  set  to  zero  and  can  be  determined  through  a 
set  of  extended  non-linear  equations  (Marchisio  & Fox  2003). 


3.  Numerics 

The  numerical  implementation  for  reactive  flow  simulations  consists  of  two  parts  - 
the  flow  solver  and  the  scalar  handler.  Here,  we  use  a LES  flow  solver  along  with  three 
different  implementations  of  the  scalar  transport  equation.  The  first  method  solves  for 
the  mean  mixture  fraction  and  reaction  progress  variable  along  with  the  variance  of 
the  mixture  fraction  using  Eulerian  transport  equations  for  these  quantities.  The  second 
method  uses  a Lagrangian  Monte-Carlo  scheme  to  solve  the  FDF  of  the  mixture-fraction 
and  reaction  progress  variable.  The  third  implementation  uses  the  DQMOM  based  2- 
environment  model  to  obtain  the  first  and  second  moments  of  the  scalars.  In  all  these 
cases,  the  density  field  is  obtained  from  the  scalar  handler  and  fed  to  the  flow  solver  which 
corrects  the  flow  according  to  reaction.  To  first  establish  the  viability  of  the  DQMOM 
model,  we  only  consider  reactions  with  no  heat  release  in  constant  density  incompressible 
flows  and  thus  do  not  use  the  feedback  loop.  The  flow  solver  provides  the  one-way  transfer 
of  velocity  and  turbulence  fields  to  the  scalar  handler.  The  following  subsections  briefly 
explain  the  numerical  implementations  of  each  of  the  components. 

3.1.  LES  solver 

The  second  order  LES  scheme  uses  an  energy  conserving  formulation  for  the  momentum 
equations  (Pierce  2001).  The  eddy  viscosity  and  diffusivity  are  computed  using  dynamic 
Smagorinsky  model  (Moin  et  al.  1991).  The  LES  scheme  also  solves  for  the  scalar  trans- 
port equation  using  a semi-implicit  scheme.  In  the  current  work,  scalar  equations  are 
solved  for  the  mixture  fraction  Z and  the  reaction  progress  variable  Y.  The  numerical 
implementation  uses  an  upwind  based  QUICK  scheme  that  is  designed  to  reduce  numer- 
ical oscillations  but  also  leads  to  some  amount  of  numerical  dissipation. 

For  the  sake  of  comparison,  the  mixture  fraction  variance  transport  equation  is  also 
simulated  using  production  and  dissipation  terms  consistent  with  the  FDF  formulation. 
Since  sub-grid  variance  is  a very  small  quantity  in  LES  simulations,  numerical  dissipation 
can  further  increase  the  errors  in  the  computation.  To  overcome  this  problem,  a ~Z  — Z2 
system  is  solved  rather  than  the  variance  transport  equation.  The  scalar  equations  are 
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obtained  by  integrating  the  FDF  equation  over  mixture  fraction  space  after  multiplying 
by  Z and  Z 2 . 

^ + VpUZ  = v(r  + rT)vz  (3.1) 

^|^  + VpUZ2  = v(r  + rT)vz5-fim(z2-zz),  (3.2) 

with  Qm  defined  as  before  providing  the  time  scale  for  scalar  dissipation.  This  formulation 
ensures  that  in  the  limit  of  zero  dissipation,  the  sub-grid  variance  is  identical  to  the 
analytical  solution: 

Zd2  ^Z2 -Z  = Z (l -Z) . (3.3) 

The  inlet  and  boundary  conditions  for  Z2  are  identical  to  the  filtered  mixture  frac- 
tion equation.  This  ensures  that  the  sub-grid  variance  is  zero  at  both  the  inlet  and  any 
boundaries  in  the  domain.  The  chemical  source  term  for  the  progress-variable  equation 
is  modeled  based  on  the  filtered  means  of  the  scalars  and  neglecting  any  sub-grid  fluc- 
tuations. This  “laminar”  assumption  leads  to  the  following  transport  equation  for  the 
reactive  scalar: 

^ + VpUy  = v (r  + rT)  vy  + pS(Y).  (3.4) 

The  LES  solution  of  the  mixture  fraction  mean  and  variance  equation  provides  an  inde- 
pendent way  of  checking  the  transported  PDF  and  DQMOM  methods. 

3.2.  Transported,  PDF  solver 

In  the  Lagrangian  implementation,  the  FDF  transport  equation  is  solved  using  an  en- 
semble of  notional  particles.  The  hybrid  composition-PDF  scheme  used  here  uses  the 
flow  fields  from  the  LES  solver  (that  was  described  in  the  previous  section)  to  evolve 
stochastic  particles  that  are  evenly  distributed  in  the  entire  computational  domain.  The 
evolution  equations  for  the  particles  are  obtained  from  Eq.  2.3  by  using  techniques  sim- 
ilar to  the  RANS  based  hybrid  method  (Pope  1985;  Colucci  et  al.  1998).  The  SDE’s  are 
solved  for  a system  of  particles  which  represents  the  FDF  in  composition  space.  The 
particles  move  in  physical  space  according  to: 

dx*  = (u  + V(r  4-  rt))  dt  + ^2(r+_FrW,  (3.5) 

where  x*  is  the  instantaneous  position  of  the  particle  and  dW  represents  the  Wiener 
diffusion  term.  The  motion  in  composition  space  is  due  to  mixing  and  reaction. 

d<t> * = [Sim  (0  - 0*)  + S(0)]  dt,  (3.6) 

where  0*  is  the  notional  particle  composition.  Here  the  composition  vector  0 = [Z,  F] 
evolves  using  reaction  source  terms  5(0)  = [0,  S(Z,Y)}.  The  use  of  the  IEM  model  with 
the  appropriate  dissipation  rate  ensures  that  the  variance  dissipation  is  consistent  with 
the  LES  solver.  Since  typical  LES  simulation  use  millions  of  computational  cells,  even  a 
particle  number  density  of  10-20  particles  per  cell  will  be  computationally  expensive.  To 
ensure  tractability,  several  novel  computational  algorithms  are  implemented.  A pointer- 
based  storage  structure  is  used  for  handling  the  particles.  Conventional  algorithms  use 
sorting  procedures  to  keep  the  particles  in  order  but  are  not  feasible  for  such  large 
numbers.  Here  an  integer  numbering  procedure  is  used  that  uses  pointer  based  linking  of 
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particles  to  cells.  Cells  store  only  the  index  of  the  particles  in  their  domain.  Ownership 
of  particles  is  transferred  to  the  destination  cell  when  particles  move  across  faces.  Such 
integer  operations  are  inexpensive,  even  though  they  tend  to  loose  some  efficiency  with 
iterations  as  the  particles  occupy  discontinuous  locations  in  memory.  Particle  motion  is 
handled  using  an  element-to-element  tracking  procedure  (Subramaniam  & Haworth  2000; 
Raman  et  al.  2003  a).  This  algorithm  splits  the  motion  into  a sequence  of  sub-steps  that 
move  individual  particles  from  face-to-face  with  intermediate  re-interpolation  of  particle 
velocity.  Reflective  boundaries  can  then  be  handled  directly.  In  all  the  simulations  carried 
out  here,  a nominal  number  density  of  30  particles  per  cell  is  used.  Particle  splitting  and 
merging  is  utilized  to  minimize  the  fluctuations  in  the  number  density.  The  details  of  the 
implementation  and  some  preliminary  results  can  be  found  in  Raman  et  al.  (2003 b). 


3.3.  DQMOM  implementation 

The  two-environment  version  of  the  model  solves  for  the  scalars  w\  (u>2  = 1 — iuj),  Gn, 
G12,  G21,  G22.  The  working  set  of  equations  can  be  written  as 

^p-  + VpuWl  = V(r  + IY)Vu;i.  (3.7) 


dpGai  — 


dt 


+ VpuGai  — V (r  + Tx)  VGaj  + Sa 


(3.8) 


It  can  be  seen  that  the  evolution  equation  for  the  environment  (Eq.  3.7)  is  identical 
to  the  mixture-fraction  equation  with  no  source  term.  The  weighted  scalar  equations 
(Eq.  3.8)  contain  source  terms  that  can  be  decomposed  into  mixing,  reaction  and  correc- 
tion components.  These  terms  can  be  specified  by  first  determining  the  location  of  the 
delta-functions  in  composition  space  in  a given  computational  cell. 
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The  mixing  source  term  takes  the  form: 

■521  = -S'™!  = Clm  (waiGa 2 — wa2 Gol)  . 

The  correction  term  is  obtained  using  the  location  of  the  delta  peaks  as 
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The  reaction  source  term  is  computed  based  on  the  composition  vector  4>i  in  each  envi- 
ronment as 

S'*  = Safa)-  (3.12) 

The  source  term  for  the  scalar  equation  Gai  is  then  given  by 


Sai  = S™  + Scai+WiSrai,  (3.13) 

where  the  reaction  source  term  is  multiplied  by  the  weight  in  order  to  be  consistent  with 
the  formulation. 

The  main  implementation  issue  comes  for  the  correction  term.  The  correction  term, 
S°ai,  compensates  for  the  excess  moment  source  term  arising  from  the  finite-peak  repre- 
sentation of  the  FDF.  As  seen  above,  this  source  term  involves  the  inverse  of  the  separa- 
tion distance  of  the  delta-peaks  composition  space.  In  regions  of  near-complete  mixing, 


Figure  2.  Instantaneous  mixture  fraction  (left)  and  progress  variable  (right)  contours 
simulated  using  the  FDF  scheme. 


this  term  will  approach  zero  leading  to  very  large  correction  terms.  Several  alternatives 
have  been  suggested  (Wang  & Fox  2003)  but  in  this  work,  an  ad-hoc  limit  is  used  such 
that  for  all  points  where  the  difference  <f>ai  — (j>a 2 is  less  than  e,  this  term  is  set  to  e.  Such 
an  implementation  using  e = 10~4  was  found  to  be  stable  for  all  the  cases  studied. 

Another  important  aspect  to  note  is  that  the  scalar  solver  does  not  maintain  the 
bounds  on  the  scalars.  This  might  lead  to  peak  locations  outside  the  accessible  range  for 
the  scalars.  It  was  found  that  any  attempt  to  limit  these  values  led  to  a sharp  decrease  in 
the  estimate  of  the  variance.  To  counter  this  problem,  for  all  grid  cells  that  contain  peaks 
outside  the  natural  limits,  the  source  terms  were  set  to  zero.  The  presence  of  reaction 
source  terms  with  stiff  kinetics  was  found  to  adversely  affect  this  occurrence.  Smaller 
time-steps  that  are  determined  based  on  reaction-time  scale  were  found  to  alleviate  this 
problem. 


4.  Simulations 

The  shear  flow  geometry  of  Mungal  & Dimotakis  (1984)  is  used  to  test  the  new  scheme. 
The  configuration  consists  of  a planar  shear  layer  formed  by  two  streams  entering  at 
8.8  m/s  and  22  m/s  velocity.  Though  the  experiment  involves  a low  heat-release  fast 
chemistry,  in  this  work  we  have  not  implemented  this  mechanism.  Since  the  purpose  is  to 
compare  different  reaction  models,  a simple  first  order  mechanism  of  the  type  A+B  — > P 
is  used  for  modeling  reaction.  The  rate  expression  is  simplified  using  a mixture  fraction- 
progress  variable  approach.  Three  different  rate  expressions  that  commonly  occur  in 
reacting  flows  are  tested.  The  transported  PDF  scheme  is  also  simulated  for  the  same 
flow  conditions  and  the  results  are  compared  with  the  DQMOM  and  LES  simulations. 
A 256  x 128  grid  spanning  80 D in  the  axial  direction  and  40 D in  the  cross-stream 
direction  is  used.  Here  D is  set  such  that  the  observation  point  in  the  experiment  is 
around  50 D.  The  inlet  velocity  profiles  are  assumed  to  be  laminar  with  flat  profiles.  A 
development  region  corresponding  to  10 D extends  into  the  domain  where  the  two  streams 
are  separated  by  a splitter  plate. 
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In  all  the  simulations,  the  rate  expression  for  the  progress  variable  is  of  the  form 


rY=K 


h~Y){r- 


Z 


Jst 


(4.1) 


The  rate  constant  K is  varied  to  study  the  effect  of  non-linearity  on  DQMOM  predic- 
tions. In  the  first  case,  the  rate  constant  is  set  to  a value  of  K = 2.  Figure  2 shows 
the  instantaneous  mixture  fraction  and  progress  variable  contours  near  the  center  of  the 
domain.  It  shows  the  vortex  like  structure  common  to  shear  flows  and  the  presence  of 
highly  mixed  reactants  at  the  center  of  such  vortices.  The  peak  in  source  term  and  mean 
progress  variable  are  observed  near  the  centerline.  The  vortices  were  found  to  stretch 
to  a maximum  of  10D  in  the  cross-stream  direction.  It  was  found  that  the  DQMOM 
model  exhibits  similar  high  reaction-rate  zones.  However,  the  LES  simulations  show  a 
thick  reaction  zone  with  maximum  allowable  reaction  rate  at  each  location.  In  order  to 
compare  the  steady-state  trends,  mean  and  variance  of  all  scalars  were  time-averaged 
for  at  least  one  flow  through  time.  Figure  3 shows  the  cross-stream  profiles  of  the  mix- 
ture fraction  and  sub-grid  variance  computed  from  all  three  schemes.  Theoretically,  the 
sub-grid  variance  obtained  from  all  these  methods  should  be  identical.  However,  the  dif- 
ferences in  the  implementation  cause  some  discrepancy.  Nevertheless,  the  time-averaged 
filtered  mixture  fraction  and  sub-grid  variance  predicted  by  all  the  schemes  are  in  good 
agreement,  thereby  validating  both  the  DQMOM  and  transported-PDF  implementation. 


The  time-averaged  mean  and  sub-grid  variance  of  the  reaction  progress-variable  ob- 
tained from  the  different  schemes  show  some  interesting  features  (Fig.  4).  The  sub-grid 
variance  is  non-zero  only  for  the  DQMOM  and  transported-PDF  schemes  and  is  set  to 
zero  for  the  LES  scheme.  In  this  context  the  LES  solver  can  be  considered  as  a one- 
environment  model  with  complete  sub-filter  mixing.  If  the  transported-PDF  scheme  is 
considered  as  a multi-environment  DQMOM  model,  the  particle  scheme  with  a nominal 
number  density  of  N represents  an  N-environment  decomposition  of  the  FDF.  The  cross 
stream  profiles  of  the  mean  show  that  the  DQMOM  method  provides  a vast  improvement 
over  the  one-environment  solution.  The  mean  profile  shows  that  in  spite  of  the  simple 
rate  expression,  second  moment  terms  cannot  be  neglected.  The  differences  between  the 
LES  and  DQMOM  models  are  highest  in  the  initial  section  where  the  effect  of  unmixed 
reactants  will  be  very  important.  Since  the  inflow  is  laminar,  the  mixing  layer  itself  does 
not  become  turbulent  until  about  X = 7.5.  However,  the  LES  solver  predicts  very  high 
reaction  rates  in  even  these  laminar  regions  where  low  mixing  should  essentially  keep 
the  reaction  rates  to  a very  low  value.  This  “early-ignition”  is  observed  in  the  profiles  at 
X = 20  where  the  mean  value  predicted  by  the  LES  solver  is  atleast  50%  higher  than 
that  predicted  by  the  DQMOM  model.  Surprisingly,  the  sub-grid  variance  profile  from 
the  DQMOM  scheme  also  shows  very  good  agreement  with  the  PDF  scheme.  This  essen- 
tially implies  that  the  third  and  higher  moments  of  the  reactive  scalar  can  be  neglected 
for  this  chemistry  scheme. 

In  the  next  case,  a more  complex  rate  expression  is  implemented.  Reaction  rates  ap- 
pearing in  combustion  have  a strong  dependence  on  temperature.  Most  source  terms 
have  an  exponential  dependence  on  local  temperature.  To  simulate  such  a condition,  the 
reaction  rate  constant  was  set  to 

K - lOOOexp  V))  . (4-2) 

where  b denotes  the  degree  of  dependence  on  temperature.  For  practical  combustion 
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Figure  3.  Comparison  of  time-averaged  mean  and  variance  of  mixture  fraction  using  a 
constant  at  different  axial  locations.  ( ) DQMOM,  ( ) LES  and  ( ) PDF. 


applications,  b is  usually  set  to  values  between  5 and  6.  However,  such  high  values  lead 
to  extinction  for  the  present  flow  configuration.  Instead,  a lower  value  of  1 is  chosen  so 
that  the  reaction  zone  can  be  anchored  near  the  splitter  plate.  This  implies  a weaker 
dependence  on  temperature  but  nevertheless  makes  the  rate  expression  non-linear.  Fig.  5 
shows  the  cross-stream  profiles  of  reaction  progress  variable.  Here  again,  it  can  be  seen 
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FIGURE  4.  Comparison  of  time-averaged  (left)  mean  and  (right)  variance  of  reaction  progress 
variable  using  a constant  rate  constant  at  (top)  X=20,  (middle)  X=30  and  (bottom)  X=50. 
( ) DQMOM,  ( ) LES  and  ( ) PDF. 

that  the  DQMOM  predictions  of  both  the  mean  and  the  variance  of  the  reactive  scalar 
are  in  good  agreement  with  the  transported-PDF  results.  The  LES  predictions  show  fast 
rates  consistent  with  the  laminar  assumption.  It  should  be  noted  that  depending  on  the 
reaction  rates,  the  complete  mixing  assumption  can  also  lead  to  quenching. 
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The  third  and  final  case  uses  a polynomial  rate  function  which  is  common  to  chemical 
engineering  assumptions.  Reduced  chemical  mechanisms  like  the  chloromethane  reactions 
(West  et  al.  1999)  may  even  involve  non-integer  moments  of  the  scalar  variable.  Here, 
the  rate  is  defined  as 


ry  = 1000 


0.0001  + Y2 
1.0 + Y 


(4.3) 


This  expression  ensures  that  the  reaction  proceeds  without  the  need  for  an  ignition  source. 
The  predicted  mean  and  variance  (Fig.  6)  indicate  good  agreement  with  the  transported- 
PDF  scheme.  This  clearly  shows  that  the  DQMOM  model  with  just  two  environments 
is  able  to  drastically  improve  the  single-environment  predictions  obtained  from  the  LES 
solver. 

In  terms  of  computational  requirements,  the  particle  based  transported-PDF  solver 
was  nearly  5 times  slower  than  the  DQMOM  scheme  even  for  such  simple  geometries. 
For  complex  configurations,  the  memory  requirements  of  a large  ensemble  of  particles 
can  further  slow  down  the  simulation.  However,  the  DQMOM  scheme  it  not  without 
limitations.  In  particular,  the  time-scale  limitations  of  a stiff-chemistry  source  term  can 
reduce  the  time-step  used  in  the  DQMOM  model  while  the  robust  chemistry  solvers 
(e.g  ISAT  (Pope  1997))  can  increase  the  computational  speed  of  the  transported-PDF 
model.  In  addition,  detailed  chemistry  mechanisms  will  require  a large  number  of  scalar 
transport  equations  and  may  eventually  diminish  the  advantages  over  the  particle  scheme. 
In  spite  of  these  observations,  present  study  shows  that  the  DQMOM  scheme  is  a viable 
altrernative  to  the  Monte-Carlo  based  transported-PDF  model  and  needs  to  be  further 
explored  for  multi-species  systems. 


5.  Conclusions  and  Future  Work 

A sub-filter  model  for  reactive  flows,  namely  the  DQMOM  model,  was  formulated  for 
LES  using  the  filtered  mass  density  function.  Transport  equations  required  to  determine 
the  location  and  size  of  the  delta-peaks  were  then  formulated  for  a 2-peak  decomposition 
of  the  FDF.  The  DQMOM  scheme  was  implemented  in  an  existing  structured-grid  LES 
solver.  Simulations  of  scalar  shear  layer  using  an  experimental  configuration  showed  that 
the  first  and  second  moments  of  both  reactive  and  inert  scalars  are  in  good  agreement 
with  a conventional  Lagrangian  scheme  that  evolves  the  same  FDF.  Comparisons  with 
LES  simulations  performed  using  laminar  chemistry  assumption  for  the  reactive  scalar 
show  that  the  new  method  provides  vast  improvements  at  minimal  computational  cost. 

Currently,  the  DQMOM  model  is  being  implemented  for  use  with  the  progress  vari- 
able/mixture fraction  model  of  Pierce  (2001).  Comparisons  with  experimental  results  and 
LES  simulations  using  a single-environment  for  the  progress- variable  are  planned.  Future 
studies  will  aim  at  understanding  the  effect  of  increase  in  environments  on  predictions. 
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